qui log using SOEP-3-Estimation.log, text replace
/* ------------------------------------------------------------------------ ** 
    C. Wunder, A. Wiencierz, J. Schwarze, H. K�chenhoff:    
    Well-being over the life span: semiparametric evidence from British 
    and German longitudinal data
    
    ESTIMATION OF SEMIPARAMETRIC REGRESSIONS

    Data organization:  person-year observations
    Data source:        German Socio-Economic Panel Study (SOEP), years 1986-2007.
** ------------------------------------------------------------------------ */ 

/*  --------( load data )-------------------------------------------------- */
version 10
matrix  drop _all
use     "SOEP.dta", clear
rename  age x
rename  gls_life_sat gls_y

/*  --------( set knots )-------------------------------------------------- */
mkspline    spline 16 = x, marginal displayknots pctile
global      N_knots = r(N_knots)            
global      N_spl = r(N_knots)+1            
matrix      knots = r(knots)'

/*  --------( define macros & prepare data )------------------------------ */
global  xvar ///
        female      disabled    hospital    educ        ln_netto    ln_hhsize ///
        german      fulltime    parttime    unempl      single      divorced  ///
        widowed     west        d_year4   	attr_in_1   attr_in_2   attr_in_3 ///
        d_year5     d_year6     d_year8     d_year9     d_year11    d_year12  ///
        d_year13    d_year14    d_year15    d_year16    d_year17    d_year18  ///
        d_year19    d_year20    d_year21    d_year22    d_year23    d_year24  ///
        one

tokenize ${xvar}

global  gls_xvar ///
        gls_`1'  gls_`2'  gls_`3'  gls_`4'  gls_`5'  gls_`6'  gls_`7'  gls_`8'  ///
        gls_`9'  gls_`10' gls_`11' gls_`12' gls_`13' gls_`14' gls_`15' gls_`16' ///
        gls_`17' gls_`18' gls_`19' gls_`20' gls_`21' gls_`22' gls_`23' gls_`24' ///
        gls_`25' gls_`26' gls_`27' gls_`28' gls_`29' gls_`30' gls_`31' gls_`32' ///
        gls_`33' gls_`34' gls_`35' gls_`36' gls_`37'

global copy_gls_xvar ///
        _gls_`1'  _gls_`2'  _gls_`3'  _gls_`4'  _gls_`5'  _gls_`6'  _gls_`7'  _gls_`8'  ///
        _gls_`9'  _gls_`10' _gls_`11' _gls_`12' _gls_`13' _gls_`14' _gls_`15' _gls_`16' ///
        _gls_`17' _gls_`18' _gls_`19' _gls_`20' _gls_`21' _gls_`22' _gls_`23' _gls_`24' ///
        _gls_`25' _gls_`26' _gls_`27' _gls_`28' _gls_`29' _gls_`30' _gls_`31' _gls_`32' ///
        _gls_`33' _gls_`34' _gls_`35' _gls_`36' _gls_`37'

global  TPF ///
        spline2     spline3     spline4     spline5     spline6     ///
        spline7     spline8     spline9     spline10    spline11    ///
        spline12    spline13    spline14    spline15    spline16    
foreach var of global TPF {
    qui replace `var' = `var'^3
}
gen spline1_sq = spline1^2
gen spline1_cu = spline1^3
                               
global spline_vars spline1 spline1_sq spline1_cu ${TPF} 
foreach var of global spline_vars{
    qui bysort persnr: egen im_`var' = mean(`var')
    qui gen gls_`var' = `var' - theta_i*im_`var'
}
drop im_*
foreach var of global spline_vars{
    qui gen _gls_`var' = gls_`var'
}

/*  --------( estimation of semiparametric model eq. in 4 )----------------- */
xtmixed gls_y gls_spline1 gls_spline1_sq gls_spline1_cu ${gls_xvar}, noconst /// 
        || _all: gls_spline2-gls_spline$N_spl, noconstant cov(identity) ///
        , technique(bfgs) emdots
est store model1b
global N_fe = e(k_f)

/*  --------( prediction: smooth function )-------------------------------- */
predict reff*, reffects level(_all)     
global all_vars ${xvar} ${spline_vars}
foreach var of global all_vars{     
    qui replace gls_`var' = `var'   
}
replace gls_one = 1
                       
local i = 1                      
local j = `i' + 1
while `i' < $N_spl {
gen re_product`i' = reff`i'*spline`j'
local i = `i' + 1
local j = `i' + 1
}
egen blups_rcoef = rowtotal(re_product1-re_product$N_knots)
drop re_product*
adjust ${gls_xvar}, generate(yfixed)    
gen EBLUP = blups_rcoef + yfixed

/*  --------( 95% Bonferroni-corrected CIs )------------------------------- */
sort x
matrix accum CTC = _gls_spline1 _gls_spline1_sq _gls_spline1_cu     /// 
    ${copy_gls_xvar} _gls_spline2-_gls_spline$N_spl, noconstant 
global N_fe_spl = $N_fe + $N_knots      
matrix nullcol = J($N_knots,$N_fe,0)
matrix nullrow = J($N_fe,$N_fe_spl,0)
matrix D = (nullrow \ nullcol,I($N_knots))
matrix estimates = e(b)             
matrix est_sig_e = estimates[1,"lnsig_e:_cons"]
matrix est_sig_u = estimates[1,"lns1_1_1:_cons"]
scalar sigma_e = exp(est_sig_e[1,1])
scalar sigma_u = exp(est_sig_u[1,1])
scalar lambda = (sigma_e^2/sigma_u^2) 
matrix Cov = sigma_e^2*inv(CTC + lambda*D)  
collapse (mean) EBLUP spline1 spline1_sq spline1_cu spline2-spline$N_spl ${xvar}, by(x)
mkmat spline1 spline1_sq spline1_cu ${xvar} spline2-spline$N_spl, matrix(C_x)
matrix Var = C_x*Cov*C_x'
matrix Var_f = vecdiag(Var)
matrix Var_f = Var_f'
svmat Var_f, names(Var_f)
gen se_f = sqrt(Var_f)
scalar z = invnormal(1-0.025/83)
gen ci_upper = EBLUP + z*se_f
gen ci_lower = EBLUP - z*se_f

/*  --------( plot of smooth function )------------------------------------ */
cap svmat knots, names(knots)           
levelsof knots, local(myknots)
twoway (rarea ci_upper ci_lower x, bcolor(gs12)) ||     ///
     (line EBLUP x, sort connect(L) lstyle(solid)) ,    ///
     legend(off) yscale(range(3 8)) ylabel(3(1)8)       ///
     ytitle(life satisfaction) xtick(20(5)100) ytick(3(0.5)8)   ///
     xmtick(`myknots', tpos(inside) tlength(*2) tlc(black))
qui log close
exit
